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ABSTRACT 

High resolution, non-oscillatory, central difference (NOCD) numerical schemes are introduced 
as alternatives to more traditional artificial viscosity (AV) and Godunov methods for solving the 
fully general relativistic hydrodynamics equations. These new approaches provide the advantages 
of Godunov methods in capturing ultra-relativistic flows without the cost and complication of 
Riemann solvers, and the advantages of AV methods in their speed, ease of implementation, and 
general applicability without explicitly using artificial viscosity for shock capturing. Shock tube, 
wall shock, and dust accretion tests, all with adiabatic equations of state, are presented and 
compared against equivalent solutions from both AV and Godunov based codes. In the process 
we address the accuracy of time-explicit NOCD and AV methods over a wide range of Lorentz 
factors. 

Subject headings: gravitation — hydrodynamics — methods: numerical — relativity 

1. Introduction 

The earliest attempts at simulating relativistic flows in the presence of strong gravitational fields are 
attributed to May and White (1966, 1967) who investigated gravitational collapse in a one dimensional 
Lagrangian code using artificial viscosity (AV) methods (VonNeumann & Richtmyer 1950) to capture 
shock waves. Wilson (1972, 1979) subsequently introduced an alternative Eulerian coordinate approach 
in multi-dimensional calculations, using traditional finite difference upwind methods and artificial viscosity 
for shock capturing. Since these earliest studies, AV methods have continued to be developed in their 
popularity and applied to a variety of problems due largely to their general robustness (Hawley, Smarr, & 
Wilson 1984a, b; Centrella & Wilson 1984; Anninos 1998). These methods are also computationally cheap, 
easy to implement, and easily adaptable to multi-physics applications. However, it has been demonstrated 
that problems involving high Lorentz factors (greater than a few) are particularly sensitive to different 
implementations of the viscosity terms, and can result in large numerical errors if solved using time explicit 
methods (Norman & Winkler 1986). 

Significant progress has been made in recent years to take advantage of the conservational form of the 
hydrodynamics system of equations to apply Godunov-type methods and approximate Riemann solvers to 
simulate ultra-relativistic flows (Eulderink & Mellema 1995; Banyuls et al. 1997; Font et al. 2000). 
Although Godunov-based schemes are accepted as more accurate alternatives to AV methods, especially in 
the limit of high Lorentz factors, they are not infallible and should generally be used with caution. They may 
produce unexpected results in certain cases that can be overcome only with specialized fixes or by adding 
additional dissipation. A few known examples include the admittance of expansion shocks, negative internal 
energies in kinematically dominated flows, 'carbuncle' effect in high Mach number bow shocks, kinked Mach 
stems, and odd/even decoupling in mesh-aligned shocks (Quirk 1994). Godunov methods, whether they 
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solve the Riemann problem exactly or approximately, are also computationally much more expensive than 
their simpler AV counterparts, and more difficult to extend the system of equations to include additional 
physics. 

Hence we have undertaken this current study to explore an alternative approach of using high resolution, 
non-oscillatory, central difference (NOCD) methods (Jiang et al. 1998; Jiang & Tadmor 1998) to solve 
the relativistic hydrodynamics equations. These new schemes combine the speed, efficiency, and flexibility 
of AV methods with the advantages of the potentially more accurate conservative formulation approach of 
Godunov methods, but without the cost and complication of Riemann solvers and flux splitting. 

The NOCD methods are implemented as part of a new code we developed called Cosmos, and designed 
for fully general relativistic problems. Cosmos is a collection of massively parallel, multi-dimensional, multi- 
physics solvers applicable to both Newtonian and general relativistic systems, and currently includes five dif- 
ferent computational fluid dynamics (CFD) methods, equilibrium and non-equilibrium primordial chemistry, 
photoionization, radiative cooling, radiation flux-limited diffusion, radiation pressure, scalar fields, Newto- 
nian external and self gravity, arbitrary spacetime geometries, and viscous stress. The five hydrodynamics 
methods include a Godunov (TVD) solver for Newtonian flows, two artificial viscosity codes for general 
relativistic systems (differentiated by mesh or variable centering type: staggered versus zone-centered), and 
two relativistic methods based on non-oscillatory central difference schemes (differentiated also by the mesh 
type: staggered versus centered in time and space). The emphasis in the following sections is to review our 
particular implementations of the AV and NOCD methods and compare results of various shock wave and 
accretion test calculations with other published results. We also explore the accuracy of both AV and NOCD 
methods in simulating ultra-relativistic shocks over a wide range of Lorentz factors. 



2. Basic Equations 

2.1. Internal Energy Formulation 

Both of the artificial viscosity methods in Cosmos are based on an internal energy formulation of the 
perfect fluid conservation equations. The equations are derived from the 4- velocity normalization u^u^ = — 1, 
the conservation of baryon number V^(pu M ) = for the fluid rest mass density, the parallel component of 
the stress-energy conservation equation u v \7 fl T tll/ = for internal energy, the transverse component of the 
stress-energy conservation equation (g av + u a u„)V ^T* 1 " = for momentum, and an adiabatic equation of 
state (eos) for the fluid pressure P = (T — l)e, where T is the adiabatic index and e is the fluid internal 
energy density. For a perfect fluid, the stress-energy tensor is 

T"" = phu^u" + Pg^ u , (2-1) 



where 



h = l + e+ — = l + re 
P 



(2-2) 
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is the relativistic enthalpy, e is the specific internal energy, u M is the contravariant 4-velocity, and g^ is the 
4-metric. The resulting equations can be written in flux conservative form as (Wilson 1979) 

d _R + d _^n o (2-3) 
dt dx* ~ ' [ ' 

dE d(EV^) dW dQVVi) 
dSj d{SjV l ) S^S V dg^ u ^-dP 



dt + dx* 25° dxi + ^8xJ ~ °' (2_5) 

where g is the determinant of the 4-metric, W = ^J~^gu° is the relativistic boost factor, D = Wp is the 
generalized fluid density, V 1 — u l /u° is the transport velocity, Si — Wphui is the covariant momentum 
density, and E = We = W pe is the generalized internal energy density. We use the standard convention in 
which Greek (Latin) indices refer to 4(3)-space components. 

The system of equations (2-3) - (2-5) are complemented by two additional expressions for V 1 and W 
that are convenient for numerical computation. Introducing a general tensor form for artificial viscosity Qij 
(see section 3.1), and defining 

,o 



M = phW = phj~gu y) =E + D + (P + ix[Qij])W, (2-6) 

the momentum can be expressed as S^ — Mu^, and So is computed from the normalization of the four- 
velocity S^Sfi = —M 2 . The coordinate velocity then becomes V % — S l /S° with V° = 1. Also, the time 
component of the four-velocity u° can be calculated from the normalization u^u^ — vPV^S^/M = — 1, and 
used to derive the following expressions for W 

W= -yF9M = y E g^_ 

S^V» phW y ' 

The former expression (W — —^p^g M/(S f _ t V' J ')) is used in the staggered mesh AV schemes as it results in 
more accurate density and velocity jump conditions across shock fronts. The latter is more convenient for 
the zone centered NOCD methods. 



2.2. Conservative Energy Formulation 

The second class of numerical methods presented in this paper (the NOCD schemes) are based on a 
simpler conservative hyperbolic formulation of the hydrodynamics equations. In this case, the equations are 
derived directly from the conservation of stress-energy 

Vj.T"" = -±= (V~gT^) ^ + r^T"° = 0. (2-8) 

Expanding (2-8) into time and space explicit parts yields the flux conservative equations for general stress- 
energy tensors 

9(^T°") d{^/=gT iu ) 



dt dx l K ' 

with curvature source terms 

£" = -yf=g T^ 7 . (2-10) 
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Substituting the perfect fluid stress tensor (2-1) into (2-9), and including baryon conservation results in the 
following set of equations 



dD diDV*) 

~dt + dx* ~ °' 



d£_ ^ d{£V l ) ^ d[^7j (g m - g m V l ) P] 
~dt * 



dx % 



dS j d(SW l ) d[^g~ {g ij - g^V*) P] 
~dt + dx l + dx 1 



(2-11) 
(2-12) 
(2-13) 



where the variables D, V 1 , and g are the same as those defined in the internal energy formulation. However, 
now 

£ = Wphu" + g 00 P, (2-14) 
S l = Wphu 1 + V^9 9 0i P, (2-15) 

are the new expressions for energy and momenta. 

It is convenient to express £ and S t in terms of the internal energy formulation variables, especially for 
initializing data 

W 2 (D 



£ = 



^v»' +r #) + < r - 1 ^» 00 #- 



g la S a + (T - 1) sf=g- 9 



E 



and reconstructing the equation of state 

w 2 wj r + (r - l)g 0Q (V=g/W) 2 ' 

where we have explicitly assumed an adiabatic gamma-law fluid. 



r - 1 



(2-16) 
(2-17) 

(2-18) 
(2-19) 



3. Numerical Methods 

Cosmos is a multi-dimensional (1, 2 or 3D) code that uses regularly spaced Cartesian meshes for spatial 
finite differencing or finite volume discretization methods. Evolved variables are defined at the zone centers 
in the NOCD, TVD, and non-staggered AV methods. In the staggered mesh AV method, variables are 
centered either at zone faces (the velocity V J and momentum Sj vectors) or zone centers (all other scalar 
and tensor variables). Periodic, reflection, constant in time, user-specified, and flat (vanishing first derivative) 
boundary conditions are supported for all variables in the evolutions. The hydrodynamic equations in both of 
the formalisms presented in §2 are solved with time-explicit, operator split methods with second order spatial 
finite differencing. Single-step time integration and dimensional splitting is used for both AV methods. The 
NOCD schemes use a second order predictor-corrector time integration with dimensional splitting, and the 
TVD approach utilizes a third order Runge-Kutta time integration with finite volume representations for 
source updates. Since the main emphasis here is on relativistic hydrodynamics, the following discussion is 
limited to presenting details relevant for the AV and NOCD schemes: the TVD method is currently only 
Newtonian. 
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3.1. Artificial Viscosity 

The order and frequency in which various source terms and state variables are updated in the AV 
methods can affect the numerical accuracy, especially in high boost flows. The following order composing a 
complete single cycle or time-step solution has been determined to produce a reasonable compromise between 
cost and accuracy: 

• Compute time step At from (3-1) 

• Store current value of boost factor W 

• Curvature 

— compute pressure and sound speed from the ideal fluid equation of state: 
P=(T- l)E/W, C s = y/TP/{ph) 

— evaluate scalar or tensor artificial viscosity Qij 

— normalize velocity and update boost factor: 

V i = S l /S", using S„S^ = —M 2 = —(D + E + PW + tr[Qy]W) 2 to first compute 5 ; 
then construct 5 M from g^ , So, and the evolved Sj\ 

and finally use equation (2-7) to define the boost factor W = —^J~^gMj (S^V 1 ) 

— update momentum, accounting for curvature: 

Sj = S^S" g lll ,j/(2S ), using second order finite differencing of g^ v 

• Artificial viscosity 

— compute pressure 

— normalize velocity, update W 

— compute pressure and sound speed 

— evaluate artificial viscosity components Qij 

— update momentum and energy equations accounting for Q^: 
E=- Eij Qtf[Vi(W') + V j (W')]/2, and Sj = -^7iQij 

• Compression 

— compute pressure 

— normalize velocity, update W 

— compute pressure again 

— update energy, accounting for compressional heating: 
E = -PV i {WV i ) 

• Pressure gradient 

— compute pressure 

— update momentum, accounting for pressure gradients: 
Sj - -y/=gVjP 

• Transport 
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— compute pressure 

— normalize velocity, update W 

— update transport terms in all variables: 

D = -V,(W J ), E = -ViiEV 1 ), and Sj = -\I l {S J V i ) 

• Boost factor 

— compute pressure and sound speed 

— normalize velocity, final update of W 

— update energy, accounting for the variation of W in time: 

E = -[P + (E i Q%/E i Qu)]w 

• Update spacetime metric components g^ v and g^ u if time dependent 



The highly nonlinear coupling of pressure and artificial viscosity to the state and kinematic variables 
through the Lorentz factor makes the relativistic equations much more difficult to solve than their Newtonian 
versions. It is for this reason that Norman & Winkler (1986) adopted implicit methods to solve the special 
relativistic equations. It is also why we have attempted to maintain a consistent and frequent update of the 
velocity normalization, boost factor, pressure and artificial viscosity throughout the cycle. 

To enforce stable evolutions, the time step is defined for all hydro methods as the minimum causality 
constraint over the entire mesh arising from the sound speed, fluid velocity, magnitude of the artificial 
viscosity coefficient, and any other physics criteria introduced in the calculation, say from scalar fields, 
radiation transport, gravity, etc... Also, since the time steps can be nonuniform, a final constraint is added 
to prevent At from increasing by more than 20% per time step. In short, 



Ai" +1 = min 



1.2 x At n 



where the superscript n refers to the discrete time level and the maximum velocity 
zones) accounts for local sound speed, fluid velocity, and viscous diffusion 



Vrr, 



max 



( " , max (J-jJ ) , Ak q2 max(|^|), Ak q2 \ £ V\ 



min (dx 2 ) 



(3-1) 

(computed over all 



(3-2) 



The Courant factor k c is typically set to < 0.6, the viscosity strength coefficient k q2 (defined in (3-8)) is set 
to 2.0 for all the problems presented in this paper, and the sound speed is defined as 



C 2 = 

3 h dp 



TP 



r(r- \)p 



(3-3) 



ph {T -1){D/W)+TP' 

for relativistic flows, where we have explicitly used the adiabatic eos in the form P oc p v . 

The artificial viscosity is implemented in a form that mimics a general imperfect fluid stress tensor 

= phu^u v + Pg» v - 2 V a^ - + «V), (3-4) 

where -q and £ are the shear and bulk viscosity coefficients, 9 = is the expansion of fluid world lines, 
and is the trace-free spatial shear tensor. Artificial viscosity is included as a bulk scalar viscosity so the 
effective perfect fluid stress energy tensor takes the form 



T"" = (p + e + P + Q)u^u» + (P + Q)g» v , 



(3-5) 
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which is equivalent to setting P — > P + Q in the momentum and energy equations 



dS. 



+ 



dx 1 



dE diEV 1 ) 



dt 



+ 



dx 1 



+ P + 



+ p 



d(p + Qj) 



2S° dxi v a dxi 
d{WV l ) ^ 9(W) 
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dx l 



dx l 



= o, 

= 0. 



(3-6) 
(3-7) 



The scalar viscosity Qi is computed as a local quantity in a dimensionally split fashion, and active only in 
convergent flows for which ViV 1 < 



Qt = (D + E + PW)Al(ViV i ) [k q2 M{\7 iV l )(l — 4> 2 ) — k ql C s ] . 



(3-8) 



The coefficients k q 2 and k q \ control the amount of quadratic and linear (in velocity) components of viscosity 
Al is a length scale set to the zone dimension, and <f> represents a limiter bounded by zero and unity 
that can be applied to reduce the effect of artificial heating and narrow the width of shock fronts. One 
could alternatively use W~ 1 Vi(WV l ) in place of ViV 1 in (3-8), which we find to be effective at preventing 
excessively large jump errors and helps stabilize solutions in highly relativistic shock tube and wall shock 
calculations. 

A more general tensor version of artificial viscosity is also implemented for convergent flows to the form 
(Tscharnuter & Winkler 1979) 



Qij = (D + E + PW)Al [k q2 V k V k Al - k ql C s 



(3-9) 



where c is a constant defined as zero or unity depending on whether the viscosity tensor should be traceless 
or not, and 5ij is the Kronecker delta. The equations for energy and momentum with a tensor viscosity are 
similar to (3-6) and (3-7) except in the way two of the viscosity terms are computed 



dE d{EV l ) 



dt 



+ 



dx % 



+ P + 



12i Qi 
Si Qi 



8Sj | 8(SjVi) S'S^dg^ 
dt dx 1 
dW . ^diWV 1 



dt 



+ P 



dx 1 



2S° dxi 
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.dP 

dx~i 



-dQ 



dx k 



fd(WV l ) d{WVi) 



V dxi 



dx 1 



0,(3-10) 
0.(3-11) 



The scalar form of artificial viscosity (3-8) is used in all the tests presented in this paper. 

The transport step is solved in a directionally split, flux conservative manner. For example, considering 
advection of the density field along the x-axis in a staggered mesh scheme, the solution to D = —\7 X (DV X ) 
is written 



D 



n+l 



DV 



At n 
Ax 



D l+1 V l+1 - DiVi 



(3-12) 



where V^+i is the face-centered velocity between zones i and i + 1, and Di is a first order monotonic Taylor's 
approximation of Di from the upwind cell center to the advection control volume center 



1 

+ 



sign 
1 

2 ~ 



1 



2 

sign 



, Vi 



Vi 



(Ax-VjAt) 
Di-i H (Vxl^ji-i 



(3-13) 



Equation (3-13) automatically detects the upwind cell from the sign of the velocity V. Here, sign(l/2, Vi) is 
Fortran notation for ±1/2, depending on the sign of Vi. High order van Leer (1977) monotonic interpolation 
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is used to reconstruct local gradients (VxD); and prevent spurious oscillations near regions of sharp gradients 

'1 /l 



(V x D)i 



si-n \ -, AD, AA-i 



/ 2AA AA_i 
I AA + AA-i + <5 



(3-14) 



The constant S <C 1 is introduced to prevent numerical overflow, and A A = (A+i — A)/ Ax are the 
mesh aligned gradients centered on the cell faces. Similar expressions can easily be derived for zone-centered 
variables on nonstaggered meshes by face-averaging the velocities, and for face-centered variables on staggered 
meshes by shifting the spatial indices and control volumes appropriately. 



3.2. Non-oscillatory Central Difference Schemes 

Considering the simplicity of equations (2-11) - (2-13), an obvious benefit of the NOCD approach is that, 
unlike the AV approach, it is not expected to be particularly sensitive to any ordering of operator updates 
since the method basically just solves a single first order operator equation with external sources. We have 
implemented two variations of this method: the first with non-staggered spatial and temporal meshes with 
second order reconstruction, and the second with time-staggered meshes in which the variables are updated 
on a mesh shifted in time to center the solution properly to second order. A summary of the solver sequence 
for this class of methods is: 

• Compute time step At from (3-1), 

redefine At — > At/2 for the 2-step, subcycled, staggered mesh scheme 

• Curvature 

— compute pressure from the ideal fluid equation of state: 

P = (r - 1)[£^/{W 2 ) - D/W]/[T + (r - l)g m (V=g/W) 2 } 

— update energy and momentum, accounting for curvature: 

£ = S° and S j = £■?, using second order finite differencing for metric derivatives 

• Flux operator 

— compute pressure from eos 

— normalize velocity and update boost factor: 

V* = S'/S , using SpS* = -M 2 = -[{£- ^f~gg m P)^f~glW\ 2 and 

S 1 = S 1 - v^gg 0t P to first compute 5°, then the boost factor W = y/=gS°/M 

— compute pressure 

— update all variables u> = (A A 

accounting for flux-conservative gradient terms in equations (2-11) - (2-13): 
Co = -\7,[luV 1 + J—gP{g ia - g 0a V 1 )] 

— if the mesh is nonstaggered in time: 

perform interpolations to recenter variables on the original nonstaggered mesh 

, .n+l _ ( 71+1 , 71+1 \ Icy , / 71+1' ,71+1' \ la 

U j - K-l/2 + ^ + l/ 2 )/ 2 + K-l/2 - C Vl/ 2 )/ 8 

• If the mesh is staggered: 

— repeat curvature and flux steps to evolve solution from t = t™ +1 / 2 to t n+1 
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— shift array indices to realign final coordinates at t n+1 with 
initial coordinates at t n by Wij t k = Wi-ij-i^-i 

• Update spacetime metric components g^ v and g^ u if time dependent 

Two essential assumptions built into this method are that the cell-averaged solutions can be recon- 
structed as MUSCL-type piece- wise linear interpolants, and that the flux integrals are defined and evaluated 
naturally on staggered meshes. Since we adopt directional splitting for multi-dimensional problems, the basic 
discretization scheme used to solve equations (2-11) - (2-13) can be derived from a simple one-dimensional, 
first order model equation of the form 

duo df(u)) „ .„ „ „. 

at + "ir = °' (3 " 15) 

where w represents any of the density energy or momentum variables, and f(u>) is the associated flux. A 
formal solution to (3-15) can be written over a single time cycle (t n — ► t n+1 ) on a staggered mesh as 



Wi + l/2(* n+1 )=^ + l/2(t n )-^ 



i r +i i r +i 

At J tn f^W^ ~ M J tn f ^ T » d7 



(3-16) 



Introducing the notation oj'j = — Wj-i, the average of the piece- wise linearly reconstructed solutions at 
the staggered positions ujj + i/ 2 {t n ) in (3-16) is given by 

^+1/2 = ^K + +i/2 + w 7+i/2) = + + ^K' - u'j+i)' ( 3 " 17 ) 
where vf +1 / 2 rc f° r to the piecewise linearly interpolated solutions from the upwind and downwind cell centers 

u j+i/2 = - (3-18) 

w7 +1/2 = Uj + ^K+i (3-19) 

Considering that the time averaged integrals in (3-16) can be approximated using midpoint values 

^ £ f(^(r))dr ~ /K(r +1 /2)), (3-20) 

immediately suggests a two step predictor-corrector procedure to solve (3-15): the state variables are pre- 
dicted at t = t n+1 / 2 by 

then corrected on the staggered mesh by 

-; + + i /2 - \w + + - - ^ [/«i 1/2 ) - /K" +1/2 )] . ( 3 - 22 ) 

where we have also substituted (3-17) for ujj +1 / 2 {t n ) in (3-16). Equations (3-21) and (3-22) represent the 
complete single cycle solution averaged on a staggered mesh. The mesh indices can be brought back into 
alignment by setting At — > At/2, performing two time cycle updates (computing w™^ 1 ^ then ^"^j 1 ) to time 
t n+1 = t n + At, and re-center the solution on the original zone positions by shifting the array indices as 
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As an alternative to mesh staggering, the solution after applying the corrector step can be reconstructed 
directly back to the nonstaggered cell centers by a second order piece-wise linear extrapolation 



n+1 _ 



1 



— " (, J- / \ _i_ ,n+l' , ,n+l \ 



,71+1' 



(3-23) 



to yield for the final single time-step solution 

-r 1 = jK-i+MHiJ-iKi-*) 



At 

2Aa; 



n+l/2 



)-/K_T 2 ) 



i 



,n+l 



n+1' 



g V J + l/2 J-l/2 



(3-24) 



We have found no substantial differences between the staggered and nonstaggered approaches in all the 
test calculations we have performed. Hence all subsequent results presented in this paper from this class of 
algorithms are derived with the nonstaggered mesh method using (3-21) and (3-24). 

One final important element of this method is that all gradients (of either the state variables oj'j or fluxes 
f'(u>j)) must be processed and limited for monotonicity in order to guarantee non-oscillatory behavior. This 
is accomplished with either the minmod limiter 



u'j = max 0, min ( 1 



LO j+ l 



(ujj+l -Uj). 



or the van Leer limiter 



1 + - Wj_i)/(Wj+i - U>j)\ 



( W J+1 - w j) I 



(3-25) 



(3-26) 



which satisfy the TVD constraints with appropriate Courant restrictions, although we note that steeper 
limiters can yield undesirable results especially in under-resolved high boost shock tube calculations. 



4. Code Tests 

4.1. Relativistic Shock Tube 

We begin testing the staggered AV and nonstaggered NOCD methods with one of the standard problems 
in fluid dynamics, the shock tube. This test is characterized initially by two different fluid states separated 
by a membrane. At t = the membrane is removed and the fluid evolves in such a way that five distinct 
regions appear in the flow: an undisturbed region at each end, separated by a rarefaction wave, a contact 
discontinuity, and a shock wave. This problem only checks the hydrodynamical elements of the code, as it 
assumes a flat background metric. However, it provides a good test of the shock-capturing properties of the 
different methods since it has an exact solution (Thompson 1986) against which the numerical results can 
be compared. 

Two cases of the shock tube problem are considered first: moderate boost (W = 1.43) and high boost 
(W = 3.59) shock waves. In the moderate boost case, the initial state is specified by pl = 10, Pl = 13.3, 
and Vl = to the left of the membrane and pn = 1, Pr = 10~ 6 , and Vr — to the right. In the high boost 
case, pl = 1, Pl = 10 3 , Vl = 0, and p R = 1, P R = 10~ 2 , Vr = 0. In both cases, the fluid is assumed to 
be an ideal gas with T = 5/3, and the integration domain extends over a unit grid from x = to x = 1, 
with the membrane located at x = 0.5. The AV shock tube results presented here were run using the scalar 
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artificial viscosity with a quadratic viscosity coefficient k q i = 2.0, linear viscosity coefficient k q \ = 0.3, and 
Courant factor k c = 0.6 (0.3 for the highest boost cases). For the NOCD method we use k c = 0.3 and the 
minmod limiter which gives smoother and more robust results than the steeper limiters in simulations of 
under-resolved highly relativistic shocks. We have carried out these tests in one, two and three dimensions, 
lining up the interface membrane along the main diagonals in multi-dimensional runs. We found it necessary, 
in order to maintain stability in the highest boost cases, to impose constraints on the pressure and energy 
density at each cycle to ensure they maintained positive values. Although this wasn't necessary at velocities 
smaller than about 0.95, we nevertheless enforced these conditions in all the calculations. 

Figures 1 & 2 show spatial profiles of the moderate boost results at time t = 0.4 on a grid of 400 uniformly 
spaced zones using the AV and NOCD methods respectively. Figures 3 & 4 show the corresponding solutions 
of both AV and NOCD methods for the high boost test using a higher resolution grid with 800 zones at 
time t = 0.36. The density under-shoot (about 30%) in Figures 3 and 4 is due to a lack of sufficient spatial 
resolution, and improves significantly by increasing the number of zones. Tables 1 & 2 summarize the errors 
in the primitive variables p, P, and V for different grid resolutions and CFD methods using the L-l norm (i.e., 
Il-^( a )lli = J2i,j.k AxiAyjAz k \a™ j k — Af jk \, where afj k and A^ k are the numerical and exact solutions, 
respectively, and for ID problems the orthogonal grid spacings are set to unity). Included in Table 1 for 
comparison are the errors reported by Font et al. (2000) using Marquina's approximate Riemann solver 
(Donat & Marquina 1996). They also tested the Roe and Flux-split approximate solvers and achieved 
similar results as with Marquina's method, so we do not include those numbers here. We find the errors 
in Table 1 are quite comparable between all three methods with convergence rates just under first order as 
expected for shock capturing methods. For the high boost case in Table 2, our errors are comparable to 
those reported by Marti & Muller (1996) for the same shock tube simulation using an extended high order 
piecewise parabolic method (PPM)(Colella & Woodward 1984) with an exact Riemann solver. However, 
we note that their published errors are for the conserved quantities (generalized fluid density D, generalized 
energy density r = phW 2 — P — D, and covariant momentum density S) rather than the primitive variables 
we report. Their results are included in Table 2. We also note that the slightly larger errors in the 3D AV 
results of Table 1 are due primarily to boundary effects (particularly at the grid corners) and not to shock 
capturing differences. In fact, errors computed only along the main diagonal are about the same for the 
NOCD and AV methods. 

Table 3 shows the mean-relative errors (defined as e rcl (a) = J2i. 3 ,k \ a Zj,k ~ ^j.fel/ 12i,j.k l^ij\fcl> wnere 
again a™^ k and A™j k are the numerical and exact solutions, respectively) in the primitive variables over 
a range of boost factors using 800 zones to cover the same unit domain. The different boost factors are 
established by systematically increasing the original value of Pl over the moderate boost case. These errors 
are also displayed graphically in Figure 5, comparing the AV and NOCD methods up to the maximum 
boost (W = 5.63 corresponding to a velocity of V = 0.984) allowed at this grid resolution, which we define 
as four cells to cover the leading post-shock density plateau using the analytic solution as a guide. The 
increasing trend (with boost) in error reflects the stronger nonlinear coupling through the fluid velocity and 
the narrower and steeper leading shock plateau found in the density plots of Figures 3 and 4. Over the range 
of shock velocities we have simulated, the errors are comparable between the AV, NOCD, and Godunov 
methods. 
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4.2. Relativistic Wall Shock 

A second test presented here is the wall shock problem involving the shock heating of cold fluid hitting 
a wall at the left boundary (x = 0) of a unit grid domain. The initial data are set up to be uniform across 
the grid with adiabatic index T = 4/3, pre-shocked density p\ — 1, pre-shocked pressure P\ = 7.63 x 10~ 6 , 
and velocity V\ = —v in it = — (1 — u), where v < 1 is the infall velocity parameter. When the fluid hits the 
wall a shock forms and travels to the right, separating the pre-shocked state composed of the initial data 
and the post-shocked state with solution in the wall frame 

V S = -S™ («) 

P2 - PlWi 

P 2 = p 2 (T - l)(Wi - 1), (4-2) 



Pi = Pi 



(4-3) 



where Vs is the velocity of the shock front, and the pre-shocked energy and post-shocked velocity were both 
assumed negligible (ei = V 2 = 0). To facilitate a direct comparison between our results and the Genesis 
code of Aloy, Ibahez, & Marti (1999) (which again uses Marquina's approximate Riemann solver) all of 
the results shown in the figures and tables, unless noted otherwise, were performed on a 200 zone uniformly 
spaced mesh and ran to a final time of t — 2.0. Also, for the NOCD methods, the Courant factor is set to 
k c = 0.3, and we use the van Leer limiter for gradient calculations, which generally gives smaller errors when 
compared to the more diffusive minmod limiter (about a 30% reduction for the lower boost cases we have 
tried). For the AV methods, we use the scalar viscosity with k c — 0.6, k q \ = 0.3, and k q2 = 2.0 for all the 
runs. 

Figures 6 & 7 show spatial profiles for the case with initial velocity Vinit = 0.9 and 200 zones for the AV 
and NOCD methods, respectively. Table 4 summarizes the L-l norm errors in both methods as a function 
of grid resolution. The values given in parentheses are the contributions to the total error in the first twenty 
zones from the reflection wall at x = 0. These numbers clearly indicate a disproportionate error distribution 
from wall heating, an effect that is especially evident in the AV results, and particularly in the density curve 
where the first two data points in Figure 6 differ significantly from the true post-shock state. Excluding 
this contribution may give a more accurate assessment of each method's ability to resolve the actual shock 
profile. 

Figure 8 plots the mean-relative errors (using 200 zones) in density, which are generally greater than 
errors in either the pressure or velocity, as a function of boost factor up to about the maximum boost 
that the AV methods can be run accurately. Although we are not able to extend the AV method reliably 
(which we define by a 10% mean error threshold, and increased sensitivity to viscosity parameters) beyond 
Vinit ~ 0.97, the NOCD methods, on the other hand, are substantially more robust. In fact, as shown in 
Table 5 and Figure 9, the NOCD schemes can be run up to arbitrarily high boost factors with stable moan 
relative errors, typically less than two percent with no significant increasing trend. These errors are generally 
smaller than those quoted by Aloy, Ibahez, & Marti (1999). However, we note that the errors for the AV 
method presented in Figure 8 and Table 5 can be improved significantly by either lowering the Courant factor 
or increasing the viscosity coefficients. For example, decreasing k c from 0.6 to 0.3, or increasing k q2 from 2 
to 3 for the case Vi n u = 0.95 reduces the L-l norm in density from 0.116 to 0.048 and 0.033, respectively. 
We have also been able to run accurate wall shock tests with the AV method at higher boosts than shown 
in Table 5 by choosing different parameter combinations (e.g., k c = 0.2, k q i = 1.0 and 400 zones can evolve 
flows with Vinit — 0.99 fairly well). However, rather than adjusting parameters to achieve the best possible 
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result for each specific problem, we have opted to keep numerical parameters constant between code tests, 
boost factors, and numerical methods. 



4.3. Black Hole Accretion 

As a test of hydrodynamic flows in spacetimes with nontrivial curvature, we consider radial accretion of 
an ideal fluid onto a compact, strongly gravitating object, in this case a Schwarzschild black hole. The fluid 
will accrete onto the compact object along geodesies, thus allowing the general relativistic components of our 
codes to be tested against a well-known analytic stationary solution. Assuming a perfect fluid in isotropic 
Schwarzschild coordinates 

'"' = - (™0 V + (i + £)V + *> + • («) 



where p = \J x 2 + y 2 + z 2 is the isotropic radius, the exact solution to this problem is dependent on a single 
parameter, the gravitational binding energy (tt ). In terms of this parameter (which we set to u = — 1 in 
our tests), the solution can be written 

W = V=SS°°V, (4-5) 

D = 7Vf <«» 

E = wr-^{p 2 vpy (4_7) 

v p 1 / 

where W is the boost factor, V p is the radial infall velocity in isotropic radial coordinates, D is the gener- 
alized density in isotropic Cartesian coordinates, E is the generalized internal energy in isotropic Cartesian 
coordinates, T — 5/3 is the adiabatic index, and Cd and Ce are constants of integration which we set to 
Cd = f and Ce = Cu/100 in the simulations. 

The computational domain for this problem is constructed to be (5M) 3 (where M = 1 is the black hole 
mass) and centered along the z-axis with — 2.5M < (x, y) < 2.5M. In the z-direction the inner boundary 
zone z — z m in is defined to lie outside the event horizon at z m i n = 1.5M in isotropic coordinates to guarantee 
all boundary zones are outside the horizon, and extends to z — z max — 6.5M along the x = y = line. 
Calculations are carried out on different resolution grids, ranging from 16 3 to 64 3 to check code convergence. 
All variables are initially set to negligible values throughout the interior domain (D = 10~ 2 , E = DCe/Cd, 
W = 1, and Si = V 1 = 0), and the static analytic solutions are specified as outer boundary conditions at all 
times. Along the inner z boundary, outflow conditions are maintained by simply setting the first derivatives 
of all variables to zero at the end of each time step. Thus fluid flows onto the computational grid from all of 
the analytically-specified (inflow) boundaries, and exits from the lower z-plane closest to the black hole. All 
results presented here were generated from simulations run until steady-state was achieved at t = 50M, and 
numerical parameters are defined as in previous tests, namely k c = 0.6, k q i = 0, k q 2 = 2.0 in the AV runs, 
and k c = 0.3 in the NOCD results. Table 6 summarizes the global mean-relative errors in both methods as 
a function of grid resolution. Figures 10 & 11 show spatial profiles of density and velocity along the z-axis 
for 32 3 and 64 3 zones for the AV and NOCD methods, respectively. 

Although the numerical results in Table 6 converge to the analytic solution with grid resolution, they 
converge at a rate between first and second order due in part to the treatment of boundary conditions 
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and time discretization errors. In particular, comparing the analytic and numerical solutions, we find that 
maximum relative errors occur near the event horizon along the inner z-boundary. For the AV method, the 
maximum relative errors for density and velocity with 64 3 zones are 9.16% and 2.49%, respectively, compared 
to global mean-relative errors of 1.36% and 0.63%. For the NOCD method, the maximum relative errors are 
24.4% and 7.42%, compared to global mean-relative errors of 2.11% and 0.14%. The global errors in both 
methods, in spite of being computed on a nonsymmetric Cartesian mesh, are comparable to those reported 
by other authors. For instance, Hawley, Smarr, & Wilson (1984b) saw relative errors of 1-3% in density 
and velocity near the horizon using an artificial viscosity code on a 25 x 10 cylindrical grid. Banyuls et al. 
(1997) saw mean relative errors of 2.67% and 0.99% using a Godunov-type method on a 200 x 5 spherical 
grid. Also, decreasing the Courant factor from k c = 0.6 to 0.2 reduces the errors in both AV and NOCD 
methods by about a factor of three, consistent with first order time discretization, and increases the rate of 
spatial convergence closer to second order. 

5. Conclusion 

We have developed new artificial viscosity and non-oscillatory central difference numerical hydrodynam- 
ics schemes as integral components of the Cosmos code framework for performing fully general relativistic 
calculations of strong field flows. These methods have been discussed at length here and compared also with 
published state-of-the-art Godunov methods on their abilities to model shock tube, wall shock and black hole 
accretion problems. We find that for shock tube problems at moderate to high boost factors, with velocities 
up to V <~ 0.99 and limited only by grid resolution, internal energy formulations using artificial viscosity 
methods compare quite favorably with total energy schemes such as the NOCD methods, the Godunov 
methods using the Marquina, Roe, or Flux-split approximate Riemann solvers, and the piecewise parabolic 
method with an exact Riemann solver. However, AV methods can be somewhat sensitive to parameters (e.g., 
viscosity coefficients, Courant factor, etc.) and generally suspect at high boost factors (V > 0.95) in the 
wall shock problems we have considered here. On the other hand, NOCD methods can easily be extended to 
ultra-relativistic velocities (1 — V < 10~ n ) for the same wall shock tests, and are comparable in accuracy over 
the entire range of velocities we have simulated to the more standard but complicated Riemann solver codes. 
NOCD schemes thus provide a robust new alternative to simulating relativistic hydrodynamical flows since 
they offer the same advantages of Godunov methods in capturing ultra-relativistic flows but without the cost 
and complication of Riemann solvers or flux splitting. They also provide all the advantages of AV methods 
in their speed, ease of implementation, and general applicability (including straightforward extensions to 
more general equations of state) without explicitly using artificial viscosity for shock capturing. 

This work was performed under the auspices of the U.S. Department of Energy by University of Cali- 
fornia, Lawrence Livermore National Laboratory under Contract W-7405-Eng-48. 
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Fig. 1. — Normalized results for the moderate boost shock tube test using artificial viscosity for shock 
capturing and 400 zones to cover a unit grid. The solution is shown at time t = 0.4. 
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Fig. 3. — Results at time t = 0.36 for the high boost shock tube test using artificial viscosity and 800 zones. 
The data points in this plot have been sampled to reduce overcrowding. Only 400 points are shown. 
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Fig. 4. — As Figure 3 but with the non-oscillatory central difference scheme. 
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Fig. 5. — Mean relative errors in density for both the AV and NOCD methods as a function of boost for the 
relativistic shock tube problem. All calculations were run using 800 zones up to time t = 0.4. The maximum 
boost limit shown in this plot (corresponding to a velocity of 0.984) is set by the constraint to resolve the 
leading post-shock density plateau by at least four cells in the analytic solution. 
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Fig. 6. — Results for the relativistic wall shock test with 200-zone resolution and infall velocity V = —0.9c 
using artificial viscosity. The solution is shown at time t = 2.0. 
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Fig. 7. — As Figure 6 but for the NOCD scheme. 
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Fig. 8. — Mean relative errors in density for both the AV and NOCD methods as a function of boost for the 
relativistic wall shock problem. All calculations were run using 200 zones up to time t = 2.0. The AV results 
can be significantly improved and brought closer in alignment with the NOCD results by simply reducing 
the Courant factor or increasing the viscosity coefficients over the standard values we have chosen for all the 
tests. 
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Fig. 9. — Density plots for different infall velocities in the wall shock test using the NOCD method. The 
resolution is 200 zones and the displayed time is t = 2.0. 
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Fig. 10. — Plots of density and velocity along the z-axis for the dust accretion problem using the AV method. 
The filled squares and open circles correspond to resolutions of 32 3 and 64 3 , respectively. The solid line is 
the analytic solution. The displayed time is t = 50M. 
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Fig. 11. — As Figure 10, except for the NOCD method. 
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Table 1: L-l norm errors in density, pressure, and velocity for the moderate boost shock-tube tests. 
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Table 2: L-l norm errors in density, pressure, and velocity for the high boost shock-tube tests. Quoted errors 
for the PPM method represent the conserved quantities, not primitive variables. 
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Table 3: Mean- relative errors in the primitive variables for different boost factors in the shock-tube test using 
an 800 zone grid. 
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Table 4: L-l norm errors for the relativistic wall shock test with infall velocity V = —0.9c. The values given 
in parentheses are the contribution of the first 20 zones to the total error. Wall heating dominates and 
greatly inflates the errors in regions near the reflective boundary, especially in the AV methods. 
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Wilson a 


8.29 x 10 


-2 












Marquina b 


9.66 x 10 


-3 


9.07 x 10" 


-3 


8.03 x 10 


-3 


0.05 


AV 


1.16 x 10 


-1 


8.51 x 10" 


-2 


5.76 x 10 


-2 




NOCD 


7.69 x 10 


-3 


6.12 x 10" 


-3 


6.74 x 10 


-3 


0.03 


AV 


1.33 x 10 


-1 


9.38 x 10" 


-2 


7.38 x 10 


-2 




NOCD 


9.40 x 10 


-3 


7.25 x 10" 


-3 


1.01 x 10 


-2 


IO" 3 


NOCD 


4.43 x 10 


-3 


2.73 x 10" 


-3 


4.60 x 10 


-3 




Marquina b 


7.20 x 10 


-3 


5.80 x 10" 


-3 


1.26 x 10 


-2 


10" 5 


NOCD 


2.09 x 10 


-3 


1.01 x 10" 


-3 


1.35 x 10 


-3 




Marquina b 


7.93 x 10 


-3 


1.00 x 10" 


-3 


7.20 x 10 


-3 


io- 7 


NOCD 


6.30 x 10 


-3 


5.59 x 10" 


-3 


1.29 x 10 


-2 




Marquina b 


9.30 x 10 


-3 


6.10 x 10" 


-3 


8.56 x 10 


-3 


io- 9 


NOCD 


5.82 x 10 


-3 


5.14 x 10" 


-3 


9.97 x 10 


-3 




Marquina b 


1.03 x 10 


-2 


6.52 x 10" 


-3 


8.13 x 10 


-3 


io- 11 


NOCD 


1.12 x 10 


-3 


8.27 x 10" 


-4 


5.08 x 10 


-4 




Marquina b 


3.40 x 10 


-2 


1.41 x 10" 


-3 


3.26 x 10 


-3 



Table 5: Mean-relative errors in density, pressure, and velocity over a broad range of infall velocities (|V"| = 
1 — v) in the wall shock test using a 200 zone grid. As noted in the text, the AV errors can be reduced 
significantly and brought closer in agreement with the NOCD results by either increasing the viscosity 
strength or decreasing the Courant factor. 



"Centrella & Wilson (1984) 
6 Aloy, Ibancz, & Marti' (1999) 



- 30 - 



Grid 


Method 


erel(p) 








16 3 


AV 


4.81 x 10" 


-2 


3.08 x 10" 


-2 




NOCD 


1.01 x 10" 


-1 


1.53 x 10" 


-2 


32 3 


AV 


2.70 x 10" 


-2 


1.34 x 10" 


-2 




NOCD 


4.55 x 10" 


-2 


3.26 x 10" 


-3 


64 3 


AV 


1.36 x 10" 


-2 


6.32 x 10" 


-3 




NOCD 


2.11 x 10" 


-2 


1.44 x 10" 


-3 



Table 6: Mean-relative errors in density and velocity for the black hole accretion problem at time t — 50M, 
where M = 1 is the black hole mass. 



